#!/bin/env python3
import ithildin as ith
import numpy as np
import matplotlib.pyplot as plt
import sys
from typing import List
from scipy.interpolate import RectBivariateSpline
from pydiffmap.diffusion_map import DiffusionMap
from pydiffmap import visualization as diff_visualization
import matplotlib

# Courtemanche1998 model, 2 ruimtedimensies, 21 variabelen

data = ith.SimData.from_stem("myokit12/myokit_2")

# Vermijd randeffecten en initialiatie
# TODO: begintijd zou in principe in de log moeten staan
def snij_randen_weg(variabele,begintijd, eindtijd):
    def begin(index):
        if index == 0:
            return begintijd
        return variabele.shape[index]//8
    def eind(index):
        if index == 0:
            return eindtijd
        if variabele.shape[index] < 20:
            return variabele.shape[index]
        return (7 * variabele.shape[index])//8
    return variabele[begin(0):eind(0),begin(1):eind(1),begin(2):eind(2),begin(3):eind(3)]

def snij_randen_weg2(variables,begintijd,eindtijd):
    newdata = dict()
    for key in variables.keys():
        newdata[key] = snij_randen_weg(variables[key], begintijd,eindtijd)
    return newdata

def onderzoek(t0, t1, aantal_punten, eps, bandwidth_normalize=False, normalize=False, k=64):
    vars = snij_randen_weg2(data.vars, t0, t1)
    # We bekijken de faseruimte, tijds- en ruimtecoordinaten zijn daarin niet van belang.
    for key in vars.keys():
        vars[key] = vars[key].ravel()
    # Verlaag aantal punten om geheugengebruik en tijdsduur te beperken.
    # replace=True strikt genomen incorrect maar verlaagt geheugengebruik
    # en wegens de hoge hoeveelheid punten weinig belangrijk
    np.random.seed(1)
    keuzes = np.random.choice(vars['u'].shape[0], aantal_punten,replace=True)
    for key in data.vars.keys():
        vars[key] = vars[key][keuzes]
        if normalize==True and np.std(vars[key]) != 0: # K is constant
            vars[key] = vars[key] - np.mean(vars[key])
            vars[key] = vars[key] / np.std(vars[key])
    metric = 'euclidean'
    if normalize=='mahalanobis':
        metric = 'mahalanobis'
        del vars['K']
    # Nu de randen verwijderd zijn en we punten in een faseruimte hebben, kunnens we proberen
    # diffusion maps te gebruiken.
    dmap = DiffusionMap.from_sklearn(epsilon=eps,n_evecs=20,bandwidth_normalize=bandwidth_normalize,metric=metric,k=k)
    # Don't use K because that variable is constant
    dmap.fit(np.column_stack((vars['gateui'], vars['gatexs'], vars['gated'], vars['Ca'], vars['gatew'], vars['gateu'], vars['Na'], vars['gateoa'], vars['Caup'], vars['gatef'], vars['gateua'], vars['gateh'], vars['gatexr'], vars['gatefCa'], vars['u'], vars['gatem'], vars['gatev'], vars['gatej'], vars['gateoi'], vars['Carel'])))
    scatter_kwargs = { 'c' : dmap.dmap[:,3]}
    diff_visualization.embedding_plot(dmap,dim=3,scatter_kwargs=scatter_kwargs)
    return dmap

# Zoals ‘Embedding given by first three DCs, coloured by fourth (BOFC model, %s points)’,
# maar plot ook de punten die niet gebruikt worden om de diffusion map op te stellen.
# Hopelijk maakt dat de figuur wat duidelijker ...
def plot_veel_punten(t0, t1, aantal_punten, dmap):
    vars = snij_randen_weg2(data.vars, t0, t1)
    # copied from 'onderzoek'
    for key in vars.keys():
        vars[key] = vars[key].ravel()
    print("raveled")
    np.random.seed(1) # determinism
    keuzes = np.random.choice(vars['u'].shape[0], aantal_punten, replace=True)
    # calculate the diffusion coordinates # TODO waarom vind Emacs+Python de o-trema niet goed? 
    var_phi1234 = dmap.transform(np.column_stack((vars['gateui'], vars['gatexs'], vars['gated'], vars['Ca'], vars['gatew'], vars['gateu'], vars['Na'], vars['gateoa'], vars['Caup'], vars['gatef'], vars['gateua'], vars['gateh'], vars['gatexr'], vars['gatefCa'], vars['u'], vars['gatem'], vars['gatev'], vars['gatej'], vars['gateoi'], vars['Carel'])))
    print("transformed")
    # plot it!
    fig = plt.figure(figsize=(6,6))
    ax = fig.add_subplot(111,projection='3d')
    ax.scatter(var_phi1234[:,0],var_phi1234[:,1],var_phi1234[:,2],c=var_phi1234[:,3],cmap='viridis')
    ax.set_title("Embedding given by first three DCs, coloured by fourth (Courtemanche model, 12)") # todo engels
    ax.set_xlabel(r'$\psi_1$')
    ax.set_ylabel(r'$\psi_2$')
    ax.set_zlabel(r'$\psi_3$')
    plt.axis('tight') # ? copied from pydiffmap
    fig.savefig("Courtemanche2D 3+1 eigenvectors, extended embedding.png")
    #fig.savefig("BOFC 3+1 eigenvectors, extended embedding (uvw, s).pdf")
    plt.show()
    plt.close()

# vanaf t < 50 komt er een staart bij (???)
matplotlib.use("TkAgg")
#dmap = onderzoek(50, 151, 7000, 0.05)
#plt.scatter(range(1,21), dmap.evals) # plot eigenwaarden (om dimensie te bepalen)
#plt.show()
#plt.close()
# oei lijkt 21 dimensionaal ... misschien is normalisatie nodig aangezien niet alle
# variabelen in dezelfde schaal werken?
#dmap = onderzoek(50, 151, 7000, 0.05, bandwidth_normalize=True)
# oeps dat deed niets, misschien wordt er daar iets anders mee bedoelt
#dmap = onderzoek(50, 151, 7000, 0.05, normalize=True)
#plt.scatter(range(1,21), dmap.evals)
#plt.show()
#plt.close()
# deze keer waren er onderbrekingen, maar nog niet echt duidelijk.  Kijk eens wat bandwidth_normalize=True
# voor impact heeft
#dmap = onderzoek(50, 151, 7000, 0.05, bandwidth_normalize=True, normalize=True) # geen invloed

# mss is de epsilon verkeerd gekozen:
dmap = onderzoek(50, 151, 7000, 40, normalize=True)
plt.scatter(range(1,21), dmap.evals)
plt.show()
plt.close()
# Voor een grotere epsilon is er 2D structuur zichtbaar.
# Te doen: wat betekent die epsilon precies?

# Als uitbreiding van de normalisatie, kan de Mahalanobisafstand
# gebruikt worden i.p.v. de Euclidische
#dmap = onderzoek(50, 151, 7000, 40, normalize='mahalanobis')
#dmap = onderzoek(50, 151, 7000, 100, k=64, normalize='mahalanobis')
#????



# hogere k, theoretisch hogere preciese, en ik zie minder 1D-stukken
# --> sterkere daling op eigenwaardeplot ook!
# ‘Futuristisch ruimteschip’ <---- 1
dmap = onderzoek(40, 151, 8000, 40, normalize=True,k=3048)
plt.scatter(range(1,21), dmap.evals)
plt.show()
plt.close()

# Cilinderachtig
dmap = onderzoek(33, 151, 5000, 70, normalize=True,k=5000)
plt.scatter(range(1,21), dmap.evals)
plt.show()
plt.close()

# 32: steekt puntje uit
# 31: kort lijntje
# 29,30: extra lijntne
# 28: koordje
dmap = onderzoek(33, 151, 6000, 400, normalize=True,k=6000)
plt.scatter(range(1,11), dmap.evals[0:10])
plt.show()
plt.close()



#varnames = ['gateui', 'gatexs', 'gated', 'Ca', 'gatew', 'gateu', 'Na', 'gateoa', 'Caup', 'gatef', 'gateua', 'gateh', 'gatexr', 'gatefCa', 'u', 'gatem', 'gatev', 'gatej', 'gateoi', 'Carel']
varnames = ['Na', 'u']
for j in range(0,3):
    for i in range(0,len(varnames)):
        fig = plt.figure(figsize=(6,6))
        ax = fig.add_subplot()
        ax.scatter(dmap.dmap[:,j], dmap.data[:,i])
        ax.set_xlabel(r'$\psi_%s$' % (j + 1))
        ax.set_ylabel(varnames[i])
        fig.savefig("2Courtemanche2D %s <-> eigenvector %s.png" % (varnames[i], j + 1))
        plt.show()
        plt.close()

# phi1 bepaalt (+-, twee oplossingen zijn ok):
#  * gateui': interval
#  * gated
#  * Ca: soms interval
#  * gatew
#  * gateu (slechte conditie in andere richting, niet natuurlijk)
#  * Na: vlak
#  * gateoa: ja (+ knoop)
#  * Caup: ja (+ knoop)
#  * gatef: ja (+ kleine opblazing)
#  * gateua: ja (+ knoopje
#  * gateh: ja  (slechte c omg,  niet nat)
#  * gatexr: interval!!
#  * gatefCa: ja (+ opblazing)
#  * u: ja (+ knoop)
#  * gatem: ja (+ s-vorm?)
#  * gatev: ja maar scherpe overgang, slechte cond van gatev -> phi1)
#  * gatej: ja (+ sl)
#  * gateoi: ja (+ gecomp hoek)
#  * Carel: ja maar beetje dik
#
# phi2:
#  * gateui: interval
#  * gaexs: tang (?)
#  * gated: beetje int, ext)
#  * Ca: ja, int
#  * [...]
#
# gateui, gateexs, gated, Ca -> phi3: driehoek, vlak .... (bah)
# phi3 <-> Na: soort-van lijn
# gateui, gateexs, gated <-> phi4: sooooort van met interval, lintje

fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111,projection='3d')
ax.scatter(dmap.data[:,1],dmap.data[:,5],dmap.data[:,2])
plt.show()
plt.close()

# 6 loodrecht:
#   * 0: driehoek
#   * 1: onderbroken vlak
#   * 2: dun vlak
#   * 3: overgang vlak
#   * 4: overgang ijl vlak
#   * 5: ijl vlak met dichte randjes
#   * 7:  beetje dun, dikke rand
#   * 8: beetje rond, ijlig
#   * 9, 19: twee stukken die overlappen?
#   * 12: idem?
#   * 13: ????
#   * 16: pingpong
# 0 loodrecht 6
for i in range(0,20):
    fig = plt.figure(figsize=(6,6))
    ax = fig.add_subplot(111,projection='3d')
    print(i)
    ax.scatter(dmap.data[:,i],dmap.data[:,6],0*dmap.data[:,19],c=dmap.data[:,0],cmap='viridis')
    plt.show()
    plt.close()

# Mss eerst S-vormige dingen preprocessen?
# Wat met 16?
for i in range(0,20):
    print(i)
    plt.scatter(range(0,len(dmap.data[:,i])), np.sort(dmap.data[:,i])); plt.show(); plt.close()

# Bij 19 verschijnt het staartje van de cilinder ---> idee, verwijder varnames[i]
# (varnames[i] is bijna discreet, snelle overgang, moeilijk om mee te werken msischien?
# Ook al is diffusion maps linear, is voortransformeren nuttig?)  Carel (,?? nee verkeerde
# variabele

# 2, 6, 19 lijkt op cil
# 2, 4, 19

# phi4, Na: vlak?
'gatexs', 'gated', 'Ca', 'gatew', 'gateu', 'Na', 'gateoa', 'Caup', 'gatef', 'gateua', 'gateh', 'gatexr', 'gatefCa', 'u', 'gatem', 'gatev', 'gatej', 'gateoi', 'Carel'

#  * gate

Ca, Carel

# u <-> phi1 (
# 

    # Na, phi0 --> overal gevuld (onafhankelijk?)
# Na, 

# Voor t0=60: dun lijntje
# Voor t0=55: bril met pootje af
# Voor t0=58: slang
# Voor t0=40: poot vervormt
#dmap = onderzoek(40, 151, 5000, 40, normalize=True,k=4000)
#plt.scatter(range(1,21), dmap.evals)
#plt.show()
#plt.close()


#dmap = onderzoek(2, 60, 5000, 40, normalize=True,k=4000)
#plt.scatter(range(1,21), dmap.evals)
#plt.show()
#plt.close()


# y = ax + 


# plot_veel_punten(50, 151, 9, dmap) # TODO te traag
#del var1, var2, var3, var4 # geheugengebruik beperken

#dmap = DiffusionMap.from_sklearn(n_evecs = 4)
#map.fit(faseruimte)


# O
